Uncertainty Modeling for Robustness Analysis of Aircraft 
Control Upset Prevention and Recovery Systems 


Christine M. Belcastro and Thuan H. Khong 1 * 

NASA Langley Research Center, Hampton, Virginia, 23661 

Jong-Yeob Shin ! 

National Institute of Aerospace, Hampton, Virginia, 23666 

Harry Kwatny 5 and Bor-Chin Chang' 

Drexel University, Philadelphia, Pennsylvania, 19104 

and 

Gary J. Balas 1 ' 

University of Minnesota, Minneapolis, Minnesota, 55455 


Formal robustness analysis of aircraft control upset prevention and recovery systems 
could play an important role in their validation and ultimate certification. Such systems 
(developed for failure detection, identification, and reconfiguration, as well as upset 
recovery) need to be evaluated over broad regions of the flight envelope and under extreme 
flight conditions, and should include various sources of uncertainty. However, formulation 
of linear fractional transformation (LFT) models for representing system uncertainty can be 
very difficult for complex parameter-dependent systems. This paper describes a 
preliminary LFT modeling software tool which uses a matrix-based computational approach 
that can be directly applied to parametric uncertainty problems involving multivariate 
matrix polynomial dependencies. Several examples are presented (including an F-16 at an 
extreme flight condition, a missile model, and a generic example with numerous cross- 
product terms), and comparisons are given with other LFT modeling tools that are currently 
available. The LFT modeling method and preliminary software tool presented in this paper 
are shown to compare favorably with these methods. 


1. Introduction 

A ircraft loss-of-control accidents comprise the largest and most fatal aircraft accident category across all civil 
^transport classes, and can result from a large array of causal and contributing factors (e.g., system and 
component failures, control system impairment or damage, inclement weather, inappropriate pilot inputs, etc.) 
occurring either individually or in combination. Research into the characterization of the aircraft loss-of-control 
phenomenon as well as loss-of-control prevention and recovery system technologies is being conducted by NASA as 
part of its Aviation Safety and Security Program (AvSSP). As shown in Figure 1, loss-of-control events can involve 
flight beyond normal operating conditions. Moreover, these conditions are not well modeled in current transport 
simulations. Validation of both the mathematical models and the systems technologies for loss-of-control conditions 
is therefore highly nontrivial. 
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Figure 1. Current Transport Simulation and Loss-of-Control Accident Characteristics 
Relative to Angles of Attack and Sideslip 

Certification of loss-of-control prevention and recovery systems (including failure detection, identification, and 
reconfiguration as well as upset recovery subsystems) for aircraft will require a comprehensive validation process 
(integrating analysis, simulation, and experimental methods) to ensure the safety and reliability of these systems. 
Robustness analysis for systems with structured uncertainty could play an important role in this process. Robustness 
is a key issue in the performance of failure detection and accommodation systems. Failure detection schemes can 
experience performance difficulties (such as false alarms) due to system uncertainties. Robustness of the control 
system can mask faults and failures and make the detection problem more difficult. It is fairly common for 
integration of failure detection and accommodation systems to be problematic if they’re designed separately. 
Robustness analysis can also identify worst-case combinations of uncertainties, faults and failures for use in guided 
Monte Carlo simulation and/or experimental studies, and could provide risk mitigation for high-risk flight testing. 
Such testing will be conducted utilizing a dynamically scaled transport aircraft that has been developed at the NASA 
Langley Research Center as part of the Airborne Subscale Transport Aircraft Research (AirSTAR) Testbed (see 
References [1] and [2]). Robustness to nonlinear parameter variations over the flight envelope and at extreme flight 
conditions must also be considered. Reference [3] provides an excellent treatment of applying robustness analysis 
methods to the clearance of flight control laws, and reference [4] provides a robustness analysis framework for 
failure detection and accommodation systems. 

Analytical robust control methods, such as the structured singular value (p, see References [3] - [4]), require 
the formulation of a linear fractional transformation (LFT) model of the uncertain system, as shown in Figure 2. 



Figure 2. Block Diagram of LFT Model 

Formulation of the LFT model can be extremely difficult and time consuming, especially for aircraft problems 
involving parametric uncertainties (see References [3] - [4] and [7] - [20]). In fact, the difficulty in formulating the 
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uncertainty model in LFT form has been a key impediment to performing robustness analyses for these systems. 
This paper presents a numerical matrix-based modeling method and preliminary software tool for computing LFT 
models from a linear parameter varying (LPV) model of the system. Section 2 summarizes the matrix-based 
computational approach and a preliminary software tool that have been developed for obtaining LFT models for 
complex systems involving parametric uncertainties, as presented in Reference [4], Section 3 discusses the 
development of LPV and LFT models for an aircraft under extreme flight conditions and provides a detailed aircraft 
example as well as several additional examples. Section 4, provides a comparison to other LFT modeling tools that 
are currently available. Section 5 presents some concluding remarks. 


2. Numerical Parametric LFT Modeling Approach 


2.1 Numerical Matrix-Based LFT Modeling Method 

The LFT model to be solved in this section is depicted in Figure 3. 
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Figure 3. Block Diagram of LFT Modeling Problem 


The matrix A(5) contains the system uncertainties, and can be represented as follows for parametric uncertainties. 


A(5) - diag [5il n j, 52ln2> • • • > Smln m ] 

(la) 

m 


dim[A(5)] = n A = Z n- 
i=l 

(lb) 

8 = [5 1 ,5 2 ,---,5 m ] G R m 

(2) 

The LFT equation associated with Figure 3 is given below. 


S(S) = L(I - AP) ~ 1 AR + Q 

(3) 

=> S(5) = Sa(S) + So 

(4) 


The matrix S(5) is a compact representation of the system model. The matrix Q represents the nominal system 
model. The interconnection matrices P, R, and L are to be determined for the uncertain component of S using the 
following equation. 

S a (§) = L(I-AP) _1 AR (5) 

Note that S A (5) contains given system matrices which are functionally dependent on the parameters 5. A solution 
for equation (5) is summarized below for S A (5) formulated as a multivariate polynomial matrix function of 5. 
Flowever, it should be noted that multivariate rational functions can also be formulated and solved using this 
approach (see Ref. [11]). 

Equation (5) can be solved for multivariate polynomial problems by replacing the matrix inversion with a finite 
series expansion and a nilpotency condition, 

S A (S) = LAR + L[AP + (AP) 2 + . . . + (AP)"] AR (6) 

(AP) r+ ‘ = 0 (7) 

where r is determined by the degree of the largest nonzero term in S A (5). An expanded definition of P, R, and L 
containing matrix partitions associated with the 8il n j blocks of A given in Eqn. (2) are defined for i, j = 1, 2,..., m. 
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L — [L 1 L 2 ' ' L m 


(8a) 


where: 
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(8b) 
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n.xn ■ 
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tirnwsXn- n.xn , 

L eR r ° WS 1 , R eR 1 cols 


i i 

and each P main-diagonal block is nilpotent of index i];: 

(P. . ) n ‘ = 0, r|j < n i; i = 1, 2, ... , m 

11 

rjj = maximum degree of 8j in S A (5) 


(9) 


(10a) 

(10b) 


The block-triangular structure of P is sufficient but not necessary for nilpotency, and other special structures can 
also be found. Solution of Eqn. (6) for the matrices L, P, R and A(5) can then be reduced to solving the following 
set of equations. 


Linear 5; Terms: 


q -Degree 5; Terms: 


Crossterms: 


L i R i =S A ’ i= !>2, ... ,m 
5 i 


L (P )^ *R =S a , i = 1, 2, ... , m ; £¥ 1,2, ... , q, 
1 “ 1 (5^ 


( 11 ) 

( 12 ) 


V 1 V 1 

in / n \ ^ 


^in T 1 


Li, (I' ) ' l> 0’ . ) z -P. . (P. . ) “T R s £. £. E, (13) 

1 1 1 2 2 2 V-l'n^ 'njn„. n T A (g; ) 1 <8; J 2 ' ’ '(g; „_ ) T 


T T T T 


1 l 2 ln T 


where: 


5 =^i +^i +- + ^i 

n T 

ii = l,2, ... ,m-(n T -l) 
i 2 = ii + 1, ii +2, ... , m - (n T -2) 


i n =ii + (n -1), ... ,m 

T T 

n T = number of parameters in the crossterm < m 


Note that the Sa terms on the right-hand side of Eqns. (11) through (13) are the known constant matrix coefficients 
associated with the indicated parameter terms in Sa(5). Moreover, depending on the number of parameters and the 
degree of each appearing in S A (5), there can be literally hundreds of S A coefficient terms and coupled matrix 
equations to be solved (or more). Moreover, satisfying these equations as simultaneously as possible to take 
advantage of any common structure (and reduce the resulting model dimension) while satisfying the nilpotency 
condition of Eqn. (7) is highly nontrivial. 

2.1.1 Solution of L, R, and Main-Diagonal Blocks of P: 

A solution for this part of the problem is given in Refs. [13] and [18], but is summarized here for completeness. 
The blocks of L and R, and the main-diagonal blocks of P are solved simultaneously for each uncertain parameter 5; 
using the linear and -degree 5j terms. The solution is accomplished such that the resulting main-diagonal blocks 
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of P are nilpotent with the appropriate index of nilpotency equal to the maximum parameter degree in Sa(5). This 
solution is accomplished numerically with a matrix singular value decomposition (svd) by recognizing that this part 
of the problem is equivalent to a 1-D state-space (minimal) realization problem and by appropriately defining the 
equivalent block Hankel matrices. The solution is accomplished for each 5, parameter as shown by the following 
theorem. 


Theorem 


Consider the linear and ^-degree 5; terms of S A (5), which can be expanded as follows 


S A (5) = [S A ] 5 ; + [S A 1 5f + ... + [S A ] 5, n ‘ 

L,c Si Si 2 Si 1 


and use the constant coefficient matrices of Eqn. (14) to form the block Hankel matrices defined below 


1 „ -naiuici i 

'5; Si 2 


where: 


Hankel 


l H 1 H 2 


H 


6i 2 - 

SA 8 i T li 

•••SA 

6^. 0] 

H, 

H„ ■•• 

1 

2 

H. 


2 


H 

0 ■•• 


H 


0 


(14) 

(15) 

(16) 

(17) 


Eqn. (16) can be constructed from (15) by shifting each block row up and filling in the bottom block row with zero 
blocks. Define the matrix svd of Eqn. (15) as follows. 


S, = U s V x 

A 0g S * i 8 i 5 i 


T =(Us i £6 i 1/2 )(2S i 1/2 VS i T )=LiRi 


(18) 


where: 


rank(SAQ ) = rank(Li ) = rank(R; ) 
°i 


Then, {P u , R } is controllable and {Li, P n } is observable, and the matrices L,. R ; , and P M form an irreducible 
realization of S A (5) as defined by Eqn. (14), where: 




Li — [inrows L) J L i , Ri — Ri 


Pii=(Li) t S Al (R 1 ) t 


"cols 

0 


(19) 

( 20 ) 


and the notation (A)' designates the pseudoinverse of matrix A. The P;; matrix is nilpotent with index q,. 

A proof of this theorem is provided in Reference [18]. 

2.1.2 Solution of P Off-Diagonal Blocks: 

The P off-diagonal blocks are each solved using the appropriate crossterms of S A (S), as defined by Eqn. (13). 
The equation to be solved for each off-diagonal block of P is a generalized linear matrix equation. The general form 
of the equation is given below for computing P n j, where n = 1, 2, ..., m-1 and j = n+1, n+2, ... , m. 


A n [ ] Pnj Bj =SA 5n5 j [ ] 


(21) 
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The matrices A n ^ , Bj, and SAg n g. ^ in Eqn. (21) are comprised of known matrices as well as matrices that have 

already been computed at this point in the solution process. The detailed equations for Block Rows 1 and 2 are 
given below. 


Block Row 1: 


Ai' 1 Prj Bj = SAg 5 [1] 
1 J 


(22a) 


where: 


Block Row 2: 


where: 


T [1] 
Ai = 


Li 

LiPn 


LiPi 


n i -i 


— [ 2 ] — [21 
a 2 p 2j B, = Sa§ g. L 

2 J 


(22b) 


(23a) 


T [21 

a 2 = 


A™ 

A™Pl2 

—[ 1 ] 

Ai P 12 P 22 


— [1] T|2— 1 

A, P 12 P 22 


T [1] 

A 2 = 


L 2 

L 2 P 22 


L 2 P: 


112-1 


Bj Matrix Structure for Each Block Row: 


W 


solved. 


B = R 

J = [r 

PR 

P " R 

J 

1 L 1 

jj j 

JJ J 

contain the constant coefficient matrices a 


Ss.s. 

s 2 

■■■ S r, 


1 j 

8.8/ 
1 J 

5.5. 
1 J 


S 2 

S 2 2 

••• s 

Sa 6 .6. [11 = 

1 J 

s/s. 

1 J 

s/s/ 

1 J 

S. 2 S. 
1 J 


S ' 

S ' „ 

■■■ S 


6, ni S. 

sAs. 2 

8.^6. 


1 J 

1 J 

1 J 


(23b) 


(24) 


(25a) 


where: i = 1, 2, ... , j-1 , j = 2, 3, ... , m ( i *j) 

[i] 


Sa 5.s. 

2 j 


[2] 


O. [11 

^ a 8_S . 

_ 2 j 

Sa 8,8^S. 
1 2 J 


Sa 5,8_S. 

1 2 j 


Sa S,8^8. 

1 2 j 

Sa I 2 ! 

sa S,8„S. 

1 2 J 




(25b) 
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Sa 8,8„S. 

1 2 j 


[k] _ 


' S k 
¥2 8 j 

S k 2 ■ 

8 8 k S. 

1 2 J 

" 8 k q i ' 
¥2 5 j ' 

S 6 2 S k S. 

1 2 j 

S 2 k 2 ' 

5 1 6 2 j 

■' S 

6 2 8 k s¥ 
2 J 

s „ ' 

s „ 

S 

s t ls k s 
6 1 ®2 S j 

5, 8 0 k S. 2 

1 2 J 

n k n j 
8, l S„ K S. J 

1 2 j 


where: k = 1, 2, ... , r| 2 


(25c) 


The off-diagonal block equations for Block Row 1, Eqns. (22), (24) and (25a) for i=l, solve all pair-wise cross- 
product terms associated with SiSj. The off-diagonal block equations for Block Row 2, Eqns. (23), (24), and (25) for 
i=2, solve all pair-wise cross-product terms associated with 5 2 8j, plus all triple terms involving 8iS 2 5j. The third 
block row equations (not shown, but constructed similarly) solve all pair-wise crossterms associated with S 3 Sj, plus 
all triples associated with 51838 , and 5 2 8 3 8j, plus all quadruple terms associated with 8i5 2 S 3 Sj. Thus, the equations 
for the i th block row solves the pair-wise crossterms associated with SjSj, plus all combinations of crossterms 
involving 81 , 5 2 , ..., 5j and 5j. Note that the solutions include all n th -degree terms. This is accomplished by the 
main-diagonal blocks of P (P n ) raised to various powers up to q, - 1, as defined by Eqn. (10b), appearing in the A 
and B matrices of Eqn. (21). The general expressions for Eqn. (21) are given in Ref. [18], but details of an approach 
for solving them are given below. 

Eqn. (21) is a generalized linear matrix equation of the form 

AXB = C (26) 

where A, B, and C are known constant matrices. Thus, solution of the off-diagonal blocks of P can be reduced to 
solving matrix equations of the form of Eqn. (26), which requires satisfaction of the following rank conditions. 

rank[A C] = rank[ A ] , rank[B T C T ] T = rank[ B ] (27) 

Satisfaction of the column rank condition (first) and row rank condition (second) of Eqn. (27) is accomplished 
through augmentation of the dimension of the appropriate parameters, 5j, in A, which translates to an augmentation 
of the associated partitions of L, R, and P comprising the A and B matrices in Eqn. (26). This is a nontrivial task, 
because columns or rows of C cannot simply be appended to A and B. For the column rank condition, the structure 
of A becomes more complicated for successive block rows of P, for higher numbers of parameters, and for higher 
parameter degrees. For the row rank condition, the B matrix structure is fixed, but is more complicated for higher 
parameter degrees. 

Another complication to performing the augmentation is that in augmenting the underlying L, R, and P matrix 
partitions of A and B, the previously obtained solutions involving these partitions must be retained, as well as the 
nilpotency (with the correct nilpotency index) of the main-diagonal blocks of P. Moreover, the augmentation 
process must be general and implementable for any number of parameters and for any parameter degree. 

The approach taken in this paper to solve the above augmentation problem is based on utilizing the structure of 
the A and B matrices (i.e., involving successive powers of nilpotent matrices), and allows an arbitrary augmentation 
to be performed independently of the coefficient matrices in C. This enables the satisfaction of the rank conditions 
of Eqn. (27) for any problem. The basic result is given in the following Lemma, which is an extension of a Theorem 
by Flalmos in Ref. [21]. 

Lemma 

Let N e R nxn be a nilpotent matrix of index q, and M e R n ' m (n>m) be an arbitrary rank m matrix such that 
rank(N q_1 M) = m. Then, the columns of M, NM, . . . , N q_1 M are linearly independent; i.e., for n > qm: 

Rank[M, NM, ... , N q l M] = qm (28) 

where [M, NM, . . . , N q_1 M] e R nxqm . H 

This lemma directly applies to the B matrix given by Eqn. (24); its dual (obtained by taking transposes) can be 
directly applied to the A matrix for block row 1, as given by Eqn. (22b). Flowever, it can be shown that the lemma 
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can also be applied to the off-diagonal block solution for any block row by reformulating the associated A matrix to 
be in the (transposed) form of Eqn. (28). Thus, arbitrary columns can be added to the columns of Lj, and arbitrary 
rows can be added to the rows of R, during the augmentation process. Moreover, an arbitrary nilpotent 
augmentation of the correct nilpotency index can be added to the main diagonal block partitions. The augmented 
matrices become: 



R 


p 

0 " 

i 

’ p 

a 


R 

L i_ 

li 

aug 

0 

p 

11 _ 



(29) 


where Pm is obtained using the general description for a nilpotent matrix given in Ref. [22]. 

The general process is to augment Li and P,, to satisfy the column rank condition for block row i, and to 
augment R, and Py to satisfy the row rank condition for block column j. In order to retain previous solutions and to 
permit arbitrary augmentations to previously computed P,, matrices, R, (L, ) is augmented with zero rows (columns) 
when L; and P;, (R, and Py ) are augmented to satisfy the column (row) rank condition. 


2.1.3 Full P-A Model Solution: 

Once the Lj, R;, P n , and Py matrices for each parameter have been determined, the full solution is assembled. 
This is a simple matter of collecting the matrix partitions together into the full L, R, and P matrices defined in Eqn. 
(8). The A matrix is also known and given by Eqn. (2), where the number of repetitions for each parameter, n;, was 
determined in solving the L,-, R ; , Pi,, and Py matrices. 


2.2 Description of Preliminary Software Tool 

The parametric LFT modeling method described in Section 2.1 has been implemented as a preliminary software 
tool. The tool has been developed for general problems involving m parameters each raised to any maximum 
degree, plus all possible cross-product terms. The software was developed in Matlab, and requires Matlab 5 or 
above, the Control System Toolbox (for the current reduction algorithms described in Subsection 2.2.3), and the 
robust control toolbox. Subsection 2.2.1 describes the data structure and syntax for the function call. Subsection 

2.2.2 describes the LFT model construction, and Subsection 2.2.3 describes a rudimentary reduction approach that is 
currently included in the software. 


2.2.1 Data Structure and Function Call 


The main function to solve for the LFT uncertainty model is “lft model”, which has the following syntax. 

[L,P,R,Q,L_full,P_full,R_full,Q_full, delta, Snew,maxtest] = lft_model(S,tol) 

The input S is an n-dimensional cell array containing the coefficient matrices, where n is the number of 
parameters. The index of cell array S relates to the term, and its content is the coefficient matrix associated with that 
term. Since the index of a Matlab cell array does not start with zero, the index of the cell array S is defined as the 
order of the term plus one. For example, consider the following: 


S(<5) = S 0 + fi }v 


+ 


s A 


l/ + 


s A 


\x 2 y 


The highest degree of X is 2 from the term X V , and the highest degree of V is 3 from the term V . If X is 
considered as the first parameter and y the second one, the order of parameters is (x, V ) . The cell array S would 
then be entered as follows. 

S{l,l}=S 0 , S{2,1} = S Aj ,S{1,4} = S Aj and S{3,2 }=S Ai 

Note that the index of a cell is equal to the order of the corresponding parameter plus one. If the order of parameters 
is (y,x), then {l,2j = S A and so on. Zero coefficient matrices, for example, S A , S A , S A , ..., can be 


left empty, and they will be filled with zeros by the software. 
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The input “tol” is optional, and specifies the tolerance on the frobenius norm of the difference between the LFT 
model and the LPV model. The default value of “tol” is le-10. The parameter “tol” will be explained more below 
when the output “maxtesf ’ is discussed. 

The outputs L, P, R are cell arrays. L{i}, R {/ } are associated with SJ n (equation (8)), P\i, / } with 
8J H and Sjl n (equation (8b)). The output Q is a matrix; Q = S 0 . The outputs Lfull, Pfull, and R_fitll are 

matrices formed from the cells of L , P , R (equation 8), and Q full = Q. 

The output “delta” is a lxn vector, where n is the number of parameters. The i th element of “delta” is n i (see 
equation (la)). 

The output “Snew” is a cell array, which contains the cells of the input cell array S and the zero coefficient 
matrices filled out by the software. If the input cell array S has zero coefficient matrices associated with terms with 
higher degree than TJ i , the maximum degree of d i (equation (10b)), the software will delete those cells and adjust 

the index of the cell array S. If the coefficient matrices associated with one or more parameters are zero (or empty) 
matrices, the software deletes all those matrices and adjusts the dimension of cell array S. The software displays 
messages that the maximum degree of a parameter is reduced and/or a parameter is deleted. The parameters 
subsequent to the deleted parameter are moved up to replace the deleted one. The adjusted cell array S is stored in 
the output “Snew”. The input cell array S in the Matlab workspace is not altered by the software. 

The output “maxtesf ’ displays the maximum result of the comparison of the original coefficient matrices of the 
LPV model with the corresponding coefficient matrices of the LFT model and the results of the nilpotency of 
P'\i,i} using the frobenius norm. If the comparison of the LPV and LFT and/or the nilpotency of P{i J j is greater 
than “tol”, messages will be displayed to notify the user. 


2.2.2 LFT Model Construction 

This section will refer to equations in this paper and in Reference [18]. Note that the notation in this paper and 
that in Reference [18] are related as follows. P n , P l2 , P 2l , P 22 in Reference [18] are equivalent to P , R , L , 
Q in this paper. P 2X in Reference [18] is equivalent to Z ; in this paper ,P n in Reference [18] to P : - in this 

paper, Pn s . in Reference [18] to Rj in this paper. 

Let n x be the number of states, tl u be the number of inputs, and n be the number of outputs. The LFT model 
outputs are calculated in the following steps: (1.) Q, (2.) main-diagonal block P u , L j , R : associated with the 
parameter 8 j , and (3.), off-diagonal block P.. , updating L i , updating R . . These steps are described below. 


Step 1 : Calculate Q 

From equation (3), Q is the nominal system and is equal to S 0 or Q = S( l,... ,1) = >S(l) 


Step 2 : Calculate Z , i? and the main-diagonal blocks of P 


The matrices L. , , and P u (i.e., the main-diagonal blocks of P ) are computed using Matlab function 

“lft svpoly” of the LFT uncertainty modeling software. The blocks P, , L. , /? associated with each parameter d i 
are calculated simultaneously. The appropriate coefficient matrices are used to form S A and S A in 

°5i \ 

equations (15) and (16), respectively. If n. , the rank of S A , is not zero, equations (18), (19) and (20) yield the 

°8i 

blocks P u e R"< /n ‘ ; Ll sR (n ^ n ', R t e R" iX ( n * +n 8 In the Matlab workspace, P n , L / . R i are stored 
as P{i,i}, Z{/}, R{i] , respectively. If n. is zero, P{i,i], Z{/}, /?{/} are left empty. 
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Step 3 : Calculate the off-diagonal blocks of P 

The off-diagonal blocks of P are computed using Matlab function “lft mvcross” of the LFT uncertainty 
modeling software. The formation of Aj"^ , B . , and S A in equation (21) is nontrivial, and equation 3.42 in 

Reference [18] was used in the implementation. in Reference [18] is split into two matrices, which are 

translated into the notation in this paper as A^ = L s P s . The structures for L s W and P s are discussed 
in detail in Reference [18]. Forming the generalized S. [1] , S, , ..., S, ... matrix partitions (as 

~ ^ SjS j ^ ^2$ j 

shown in equations (25 a,b,c) for two parameters) for any number of parameters raised to an arbitrary order was also 
challenging. Since each cell of cell array S is associated with a term of the LPV model, allowing the number of 
parameters to change from problem to problem means that the dimension of S must change accordingly. In 
addition, since the order of the uncertain parameters change from problem to problem, the maximum index of S also 
changes. Placing the correct cell into the appropriate location in the matrix partitions of equations (25 a,b,c) was 
accomplished using the following observation in order to avoid the use of Matlab string operations. Any cell in a 
Matlab multidimensional cell array can be referred to by two methods: either by the full index of the cell or the 
equivalent one-number index of the cell. Cells in a cell array in the Matlab workspace are stacked ‘column-wise’. 
Cell \k x , . . . , k n } of cell array S {m ] m n } , (k i < m i ) can be referred to as cell \K } where K is calculated as 
follows. 

* ( i - 1 > 

k = k i + x ft -on*, 

i=2 y j = i ) 

The off-diagonal blocks of P are then calculated in the following order P n , P u , ..., P ln , P 23 , ..., P 2n , ..., 
P [n _ l)(n) . A different order of calculation of the off-diagonal blocks of P can be followed, but equation (21) and 
related equations must be revised accordingly. Equation (21), which has the form of equation (26), will yield P l} , 
i = — 1, j = i + 1 ,. . .,n if the rank condition in equation (27) is satisfied. If the column (row) rank 

condition for P tJ is not satisfied, A^ (Bj) in equation (21) is augmented to satisfy the column (row) rank 
condition. As mentioned in section 2.1, the augmentation must retain previously obtained solutions as well as the 
nilpotency of the main-diagonal blocks of P . The augmentation is described in equation (29) and in the last 
paragraph of section 2.1. In this preliminary software implementation, random matrices are generated for the 
arbitrary matrix augmentations. A more sophisticated approach will be considered for future work. 

2.2.3 One-Dimensional Model Reduction 

A modest approach to LFT model reduction is included in the preliminary software tool. Once the LFT model is 
formed, each sub-block of A is treated as a one-dimensional system and “uncontrollable” and “unobservable” modes 
are removed by performing the transformation provided in the Control System Toolbox. The one-dimensional 
system ‘A’, ‘B’, and ‘C’ matrices are formed as follows. 

A k -P{k,k}, k = \,...,n 

B k =[R{k } P{k,k + 1} P{k,k + 2} ... P{k,n}\ 

C k =[L{k\ P{\,k\ P{2,k\ \ P{k-l,k }] 

The A k , B k , C k matrices are reduced using the controllability and observability rank tests, and the reduced 
matrices are denoted as A r k , B rk , C rk . Let P rr = [rank, rank-, . . . rank n ] , where rank, is the result of the 
rank test for the i th parameter. Then the matrices L , P , R are updated as shown below. 
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P{k,k) = A rk 

L{k] — C r ^(l : « ,:) where n rows is the number of rows of a coefficient matrix 

^W = 5 ^(h 1: « C0 J where n cols is the number of columns of a coefficient matrix 

Let C b,k = C r,ki n rows + 1 : enC ^') and B l,k - B r,k (* 5 U cols + 1 : eYld ) 

P{i,k} = C Kk £ P size (1 : i - 1) + 1 : £ ^ (l : i)ll : rank, k \ i = 1,2, . . . , k - 1 
./} = % (l : rank k , £ P lze {k + 1 : j - 1) + 1 : ^P size {k + 1 : j)\ j = k + 1 , k + 2, . . . , n 


This reduction process is described in more detail in Reference [11]. Other more sophisticated reduction methods 
(see References [15]-[17]), will be incorporated into the next version of the software. 

2.2.4 Potential Software Improvements 

The software tool described in this paper is preliminary, and a number of potential improvements in its 
sophistication can be made. Some of these are considered below. 

• Numerical conditioning issues have not been fully considered. A more thorough evaluation and 
treatment of this should be performed. 

• Generating a random matrix for the augmentation was an initial simplistic approach to generating the 
arbitrary augmentation matrices for solving the off-diagonal blocks of the P matrix. A more 
sophisticated approach should be considered - perhaps one that makes use of the minimality criteria 
defined in References [15] - [17]. 

• The current software tool does not include much logic for tailoring the set of equations to be solved to 
the specific problem. Specific potential areas of improvement include the following. 

- Simplification of the equations may be possible (and preferable relative to model dimension) for 
simpler problems. 

- The software currently uses an upper block triangular structure for P - which is sufficient but not 
necessary for nilpotency. Allowing more flexibility in this structure could reduce the resulting 
model order. For example. Reference [18] includes a problem for which an LFT model of 
dimension 18 can be achieved when the P matrix is not constrained to an upper block triangular 
structure (as compared to a dimension of 19 with the constraint). 

• The structure of the equations being solved for the cross-product terms was formulated such that the 
matrix A in equation (26) absorbs most of the cross-product term complexity. It may be possible to 
reformulate the equations to balance the complexity between the A and B matrices of equation (26), 
which could possibly result in lowering the LFT model dimension. 

• A more sophisticated model reduction method, such as that presented in References [15] - [17], should 
be incorporated into the software. 

• Scaling is not currently included in the software to produce a normalized LFT. This function should be 
incorporated into the software. 


3. LFT Modeling Examples 

Several example problems are presented in this section to illustrate the development of LPV and LFT models 
using the matrix-based approach developed in this paper. Section 3.1 presents the development of these models for 
an aircraft at extreme flight conditions. Section 3.2 presents a missile model problem, and Section 3.3 provides a 
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generic example developed to be a difficult problem involving three uncertain parameters and all associated cross- 
product terms. 


3.1 Uncertainty Modeling for Extreme Aircraft Flight Conditions 

To illustrate uncertainty modeling for an extreme flight condition, LPV and LFT models were developed for an 
F-16 aircraft near a stall bifurcation. The following sub-sections describe the development of these models. Section 

3.1.1 describes the general approach for developing LPV models for extreme flight conditions, and section 3.1.2 
provides the F-16 example problem. 


3.1.1 Formulation of LPV Models Near Bifurcation Points 


An LPV model can be formulated for extreme flight conditions at (or very near) bifurcation points using the 
method proposed in Reference [24]. Consider the parameter dependent nonlinear system 


x = f(x,u,ju) 
y = g(x,u,ju) 


(30) 


Where x e R" is the state, u e R p is the control, y e R"‘ is the output, and // e R k is a parameter vector. 
Parameters may include quantities that define the desired operating condition, like speed, flight path angle, altitude 
or physical parameters such as weight, center of mass location, etc. We are interested in how the vehicle behaves 
given different parameter values. To this end we wish to construct a family of linear time-invariant models 
corresponding to different values of the parameters, // . 

A linear time-invariant model is obtained from equation (30) by linearizing at a specified equilibrium point. 
Consequently, we first need to characterize the dependence of equilibrium points on the parameters. An equilibrium 
point for the system of equations (30) is a triple (x 0 ,w 0 ,// 0 ) that satisfies the conditions f (x 0 ,u 0 ,/U 0 ) = 0 and 

g[x 0 ,u 0 ,ju 0 ) = 0 . The set of equilibrium points is a k dimensional manifold in the space R" +m * k . Formally, the 
equilibrium manifold is defined by 


6 = j(x, u, /u) e R n+m+h | / (x, u, /u) = 0, g (x, u, //) = 0] 


(31) 


The manifold will typically be quite complex and may not even be smooth. An example of a two parameter surface 
is shown in Figure 4. At any equilibrium point, it is possible to derive a linear approximation to the system (30) in 
the form 


Sx = ASx + Bu 
8y = CSx + Du 


(32) 
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The parameters A,B,C,D are obtained by evaluating the appropriate Jacobian matrices at a specified equilibrium 
point p e & . 

It is evident from Figure 4 that any attempt to express equilibrium values for x,u as functions of the parameters 
// is inherently limited by the folds and creases in the surface - i.e., the bifurcation points. An alternative approach 
can be defined in which the surface itself is parameterized so that it is possible to define functions in these new 
parameters. These functions may be globally valid, but even if local, the domain need not be constrained by 
bifurcation points. In Reference [24] it is shown how to obtain a local parameterization of & around any point 
p n e & . In this method, a parameter vector s e R k is introduced so that points p e&t cz , where &l is a 
neighborhood of p 0 in S , are defined parametrically by p = p ( s ) on a neighborhood of the origin in R k . In this 
representation, the equilibrium state, control, and parameters are all given as functions of the new parameters s . 

x 0 =x 0 (s), u 0 =u 0 (s), p 0 =p 0 (s) 


Once this parameterization of the equilibrium surface is obtained, an LPV model of the form 

Sx = A(s)<5x + .5(s)m 
Sy = C(s)Sx + D(s)u 


(33) 


is easily constructed. A software tool has been developed in Mathematica® for accomplishing this. 

3.1.2 F- 16 Example 

The LPV and LFT models for an F-16 example near a stall bifurcation are given in this section. Figures 5 and 6 
show a portion of the equilibrium surface of an F-16 near the stall condition. Details about the dynamical behavior 
of the aircraft near the bifurcation point can be found in Reference [25]. 



Figure 5. A portion of the equilibrium surface 
for an F-16 is shown. The parameters are speed 
V and flight path angle y . Only one of the 

control variables, elevator deflection 8 e , is 
shown. 


ril 



Figure 6. Another representation of the 
equilibrium surface. Slices through the surface 
at constant flight path angle are shown. From left 
to right y - -0.005,-0.0025.0,0.0025,0.005 . The 
surface clearly shows stall as speed is reduced. 
Notice that the stall speed increases with 
increasing flight path angle. 


The functions defining the equilibrium surface, x 0 (s),u 0 (s),p 0 (s) as well as the matrices A(s),.B(.s),C(.s),-D(.s) 
are obtained as polynomials in the parameters s . 


F-16 LPV Model 
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The LPV model consists of matrices A(s 1 ,s 2 ),B(s 1 ,s 2 ),C(s 1 ,s 2 ),D(s 1 ,s 2 ) . The parameters s l ,s 2 are related to 
speed, V , and flight path angle, y , through well-defined formulas, as shown below. 

V ~ fri. S Y ~ fy{. S \’ S l) 

a = 0.873652 + 0.124836 s, - 0.0607641 sf + 0.0320992 Sl 3 - 9.60292 s 2 
+ 10.371 sis 2 - 8.14459 sfs 2 - 438.22 s 2 2 + 690.549 SiS 2 2 - 19536.8 s 2 3 

Th = 15064.284 + 2256.19 s 1 - 1040.09 s, 2 + 745.335 Sl 3 - 187438 s 2 
+ 180989 sis 2 - 186906 si 2 s 2 - 7.8323xl0 6 s 2 2 + 1.56862xl0 7 SiS 2 2 - 4.4026x10 s s 2 3 

dele = 0.0827526 + 0.0717976 s, + 0.0226667 sf + 0.0239419 S; 3 - 5.52297 s 2 

- 2.89913 sis 2 - 5.00515 si 2 s 2 + 88.8865 s 2 2 + 349.239 SiS 2 2 - 8145.83 s 2 3 

0 = 0.863652 + 0.124836 si- 0.0638262 Sl 2 + 0.0381331 si 3 - 10.6029 s 2 
+ 10.8421 sjs 2 - 9.71647 s, 2 s 2 - 460.435 s 2 2 + 825.264 Sl s 2 2 - 23364.5 s 2 3 

V = 131.475 + 1.0 si+ 6.38958 si 2 - 5.1431 si 3 + 0.0 s 2 
- 954.805 SjS 2 + 1300.36 Si 2 s 2 + 35898.1 s 2 2 - 108465 Sl s 2 2 + 2.99105xl0 6 s 2 3 

y = 0.01 + 0.0 si + 0.00306217 sf - 0.00603397 Si 3 + 1.0 s 2 

- 0.471 1 1 sis 2 + 1.57187 s, 2 s 2 + 22.2155 s 2 2 - 134.715 SiS 2 2 + 3827.73 s 2 3 

Thus, each pair (spS,) corresponds to a unique (V,y) . The origin (s x ,s 2 ) = (0,0) corresponds to a flight condition 
close to the stall bifurcation point. In particular, (.s, ,.s 2 ) = (0,0) — » (V,y) = (131.45 fps, 0.01 rad) . Functions for V 
and y are plotted in Figures 7 and 8. 

V fps g- rad 



The LPV model for this example is given as: 

x = A(s) ,s 2 )x + B(s j , s 2 )u 
y = C(s t ,s 2 )x + D(s 1 ,s 2 )u 


where the parameter dependent matrices are built as follows: 

^(^,5 2 ) = AA00 + AA10s l + AA01 s 7 + AA\ l5,^ 2 + AA20s* 4 — 
and similarly for the others. 
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The LPV model has 9 states: x = {<j>,0, y/,p, q,r,V, a,j3 j , 5 controls u = { 8 a , S el ,5 er ,8 r ,Th} , and 10 outputs 
y = {x,y} , and it does capture the bifurcation behavior shown in Figures 5 and 6. The detailed F-16 LPV model is 
provided in the Appendix. 

F-16LFT Model 

The LFT can be obtained using the LFT modeling tool described in Section 2 by assigning the S matrix cell 
array as follows: 

S{1,1} =[AA00BB00], S {2, 1 } = [AA10 BB10], S{1,2} = [AA01 BB01] 

S{2,2} = [AA1 1 BB1 1], S {3,1 } = [AA20 BB20], S{1,3} = [AA02 BB02] 

S{3,2} = [AA21 BB21], S{2,3} = [AA12 BB12], S{4,1} = [AA30 BB30] 

S{1,4} = [AA03 BB03] 

Note that matrices C and D need not be included in S, since they do not contain uncertain components for this 
example. The value of tol was set to 1 X 10" 8 . The resulting LFT model dimension was 67, with 27 occurrences for 
Si and 40 occurrences for s 2 . These results can be further reduced using more sophisticated reduction methods, as 
shown in Section 4.2. 

3.2 Missile Model Problem 

The short period dynamics of a missile are considered in the example from Reference [27], with states angle of 
attack ( a ) in radians and pitch rate ( q ) in radians/sec, and control input (fin deflection) ( 8 ) in radians. The 
nonlinear equations of motion are given as follows. 

. 0.1 p a CM 2 S cos 2 (a) 

a = — 2 — + q 

( mass)U 

. 0.7 p,.C n M 2 Sd 

q = — 

where the aerodynamic coefficients are given by: 

C n = a n a 3 + b„a 2 + cf 2 + yj a + dJS 

Cm = a ,« a 3 + b ,y - cjj - ja + d m S 

U = M*ss*cos(a) 

The properties of the missile are defined as follows. 

p Q = 913.31b/ ft 2 Pressure at 20000 ft 

ss = 1031.1 ft / s Speed of sound at 20000 ft 

S = 0.44 ft 2 Reference area 

d = 0.75 ft Diameter 

mass = 13. 98 slug Mass 

M = 2 Mach number 
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The coefficients are given by: 


f = 182 .5 slug ■ ft 2 Pitch moment of inertia 


a n = 0.000 103deg“ 3 , b n = -0.00945 deg 2 , c„ = -0.1696deg"' , d„ = -0.034deg"‘ 
a m = 0.0002 15 deg" 3 , b m = -0.0195deg" 2 , c m = 0.051deg“\ d m =-0.206deg"' 


The nonlinear equations are written in quasi-LPV form as follows. 


a 

q 


0.0207 


1.232 


( mV 

M 

( 2 A 

a 


1- — 

V 3 JJ 


l 2 J 


+ b,„a - c . J 7 — M 
3 


M 2 


OmOldM 


f « 2 V 

1 -— 

V 2 J 

\232dM 1 


The nominal values of M and a are 2 and 0 radians, respectively. 

The LFT model obtained for this example has a total dimension of 10 (with 4 occurrences for a and 6 
occurrences for M). These results are compared to other software tools in Section 4.2. 


3.3 Generic Example 

In this section, an extension of a physics-based model is considered. The LPV model for this example involves: 
three uncertain parameters, one parameter with maximum degree 3 and the other two parameters with maximum 
degree 2, and all associated cross-product terms. The LPV model is given below. 
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The LFT model obtained for this example has a dimension of 45 (with 9 occurrences for x, 24 occurrences for y, 
and 12 occurrences for z). These results are compared to other software tools in the following sub-section. 


4. LFT Modeling Tool Comparison 

In this section, LFT models are obtained for the examples of Section 3 using two available software tools 
(developed by ONERA and MuSyn, Inc.) and the results are compared to those obtained using the matrix-based 
numerical software tool presented in this paper. The LFT modeling tool developed by ONERA (see Reference [23]) 
and that developed by MuSyn, Inc. (see Reference [6]) are considered, and a brief description of each is provided in 
subsections 4. 1 and 4.2, respectively. The comparison of results is presented in subsection 4.3. 

4.1 LFT Modeling Tool Developed by ONERA 

The Linear Fractional Representation (LFR) toolbox was developed by J-F Magni (see Reference [23]) based on 
an object-oriented realization technique (see Reference [28]). The description provided herein is taken from 
Reference [23]. In the toolbox, the following LFT modeling techniques have been implemented: Morton’s method 
(see Reference [8]), Horner factorization (see Reference [29]), and tree decomposition (see Reference [30]). The 
tree decomposition method is recommended in Ref. [23] as being the most efficient. The software also utilizes a 
symbolic method (based on Maple, or the symbolic tool within MATLAB) for algebraic manipulation. Morton’s 
method can represent linear parameter dependent systems as an LFT model using singular value decomposition. 
This method is applied to a polynomial parameter dependent system, including rational functional forms of the 
parameter, which results in the A matrix containing rational functions of the uncertain parameters (instead of the 
parameters themselves). The Redheffer’s star product is then used to obtain the usual form for A. Horner 
factorization concerns single variable polynomials, and its objective consists of avoiding calculation of all the 
powers of each uncertain parameter. The tree decomposition is a generalization of the idea consisting of factorizing 
parameters so that they appear a minimum number of times as is possible before proceeding to the realization. For a 
detailed description of these methods, the reader is referred to Ref.[23] and its references. For the examples 
presented in this paper, the tree decomposition method was used by invoking the “symtreed” command to construct 
the LFT models. 

The LFR toolbox also contains three different order reduction methods (see Reference [23]). The first method is 
a 1-D reduction approach, which consists of considering that each uncertain parameter plays the role of 1/s. For 
each parameter, a balanced realization is applied to reduce its size. The second approach is an n-D reduction method 
developed by Carolyn Beck et. al. (see References [15] - [17]), which consists of considering controllability and 
observability of the uncertain system. When the system matrices (A, B, C, D) have an uncontrollable or 
unobservable space, the algorithm can calculate a transformation to produce a null block on the B matrix 
(uncontrollable) or a null block on the C matrix (unobservable). Generally, the n-D approach is less conservative 
than the 1-D approach, since it treats simultaneously all the parameters of the uncertain block [17]. The third 
reduction method is the generalized Gramian approach, which consists of considering a generalized Gramian 
defined in an LMI system as follows. 


A > 0, 

M T n XM n -X+M^Mn < 0 
7 > 0 

M n YM T n - Y + M 12 Mj2 < 0 
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After using the similarity transformation, singular values less than the given tolerance are truncated and the order of 
the model is thereby reduced. A detailed description is provided in Ref. [23] and its references. 

For practical use of each reduction method, note that the n-D reduction method is computationally fast and the 
user should define a tolerance level of uncontrollability and unobservability. When the tolerance level is high (e.g., 
le-4 instead of 1 e- 12), the reduced LFT model may not capture the original dynamics. The user should therefore 
define the appropriate tolerance level for each problem. Using the Gramian approach, the user can approximate the 
system based on singular values of the Gramian (X and Y). Note, however, that the computation of the Gramian can 
have a high computational cost (calculation time) to solve the LMI optimization (see Reference [23]). 

In this paper, the n-D reduction method is applied to the LFT models obtained for the examples using the three 
modeling tools being compared (i.e., the NASA, ONERA, and MuSyn tools). The 1-D reduction approach is also 
applied to the ONERA results for comparison to the 1-D results obtained using the method presented in this paper. 


4.2 LFT Modeling Tool Developed by MuSYN, Inc. (and Included in the Matlab 7.0.4 Robust Tool) 

The following description about building uncertainty models using the Robust Control Toolbox in MATLAB 
7.0.4 is taken from Chapter 6 of Ref.[6]. In this paper, we used the Robust Control Toolbox to generate LFT 
models of the examples, using the function “ureal”. The function generates uncertain atoms for uncertain real 
parameters. Uncertain atoms are used to form uncertain matrix objects and system objects. Note that each 
uncertain atom is written in an LFT block. An LFT model of a matrix object is built up from uncertain atoms, 
depending on the sequence of operations in their construction. Note that different ways of matrix construction in 
terms of the uncertain atoms may generate different sizes of uncertainty blocks in the LFT models. 

In the Robust Control Toolbox, there are several reduction methods for LFT models. The command “simplify” 
or “AutoSimplify”, as an option to the function “ureal”, can reduce the uncertainty block size. The AutoSimplify 
parameter can be set to “off’, “basic”, or “full”. In the “off’ case, no simplification is attempted. In this paper, the 
models obtained for the “off’ case are referred to as having “no reduction”. In the “basic” case, fairly simple 
schemes to detect and eliminate non-minimal representations are used (such as removal of zero rows and columns), 
and in the “full” case, numerical-based methods similar to truncated balanced realizations are used, with a very tight 
tolerance to minimize error. The AutoSimplify property of each uncertain atom dictates the types of computations 
that are performed to generate an LFT model of the uncertain matrix or system. 

4.3 Comparison of LFT Modeling Results 

The LFT results obtained using the three software tools for the example problems presented in Section 3 are 
summarized in Table 1. It should be noted that each of the above methods produces an excellent representation of 
the given LPV model both before and after model reduction”. As indicated in the Table, the numerical matrix- 
based method and software tool presented in this paper produced an LFT model (before reduction or after 1-D 
reduction) that was comparable or lower in dimension for each of the examples than the other two tools. Models of 
substantially lower order (prior to reduction) were obtained using the NASA and ONERA tools, as compared to the 
MuSYN tool, for the F-16 and Generic examples - which involved the most complex LPV models. Using the n-D 
reduction method developed by Carolyn Beck (see References [15] - [17]) and included in the ONERA tool, the 
LFT models can be reduced to about the same low order despite what modeling tools were used. The exceptions are 
that the MUSYN “full” reduction option produced a lower-order LFT for the F-16 example (by 2 parameter 
repetitions) and a higher-order model for the Missile example (by 5 parameter repetitions). 


” The H-infinity norm was computed for each model at corner points associated with the uncertain parameters and 
over a large frequency range (0.001 - 100,000 rad/sec). The maximum deviation of the LFT and LPV models was 
on the order of lxlO" 5 . 
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Table 1. Comparison of LFT Modeling Results 


\ Example: 

F-16 

(V,y) 

Missile 
(a, M) 

Generic 
(x, y, z) 

Tool: 

n A 

n v 

n y 

n A 

n a 

n M 

n A 

n x 

n y 

n z 

NASA (“lft_model”) 











No Reduction 

86 

31 

55 

12 

4 

8 

94 

9 

64 

21 

1-D Reduction (NASA) 

67 

27 

40 

10 

4 

6 

45 

9 

24 

12 

n-D Reduction 
(ONERA) 

55 §§ 

22 

33 

9 

4 

5 

45 

9 

24 

12 

ONERA (“symtreed”) 











No Reduction 

72 

24 

48 

11 

6 

5 

110 

9 

24 

77 

1-D Reduction 
(ONERA) 

67 

23 

44 

9 

4 

5 

109 

9 

24 

76 

n-D Reduction 
(ONERA) 

56 

23 

33 

9 

4 

5 

46 

9 

24 

13 

MUSYN (“ureal”) 











No Reduction 

560 

280 

280 

16 

8 

8 

738 

318 

210 

210 

“Basic” Reduction 
(MUSYN) 

235 

106 

129 

16 

8 

8 

579 

268 

173 

138 

“Full” Reduction 
(MUSYN) 

53 

22 

31 

14 

8 

6 

45 

9 

24 

12 

n-D Reduction 
(ONERA) 

53 

22 

31 

14 

8 

6 

45 

9 

24 

12 


; Note : Reversing the order of V and y parameters resulted in a dimension of 54 (with n Y = 22 and n v = 32) after n- 
D reduction (using the ONERA tool). 
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The results reported in Table 1 represent the best results that we obtained using the MuSyn and ONERA tools. 
For the ONERA tool, the tree-decomposition method resulted in the LFT model of lowest dimension, whereas the 
“lfrs” command produced the same LFT model dimension as the MuSyn tool with the “off option. Moreover, both 
methods (i.e., the “lfrs” in the ONERA tool and “ureal” in the MuSyn tool) yield LFT models whose reduction 
depends highly on how the LPV model is represented in terms of the uncertain parameters. That is, the 1-D and n-D 
reduction results using the ONERA tool and the “full” and “basic” reduction options in the MuSyn tool appear to be 
highly dependent on the equation representation of the LPV model. To illustrate this, the Generic Example was 
equivalently implemented using the following three representations. 

LPV Representation 1 : Using a For-loop statement, the equations were written in low to high order with the 

parameter multiplications occurring in alphabetical order, as follows. 


LPV1 = zeros(3,6); 
fork= 1:3 
forj= 1:3 
for i = 1:4 

if ~isempty(S{i,j,k}) 

LPV 1 = LPV 1 + S{i,j,k}*x A (i-l)*y A (j-l)*z A (k-l); 
end 
end 
end 
end 

LPV Representation 2 : Without using a For-loop, the equations were written in ascending order with alphabetical 

parameter order, such as x,y,z, as follows. 

LPV2 = S{3,l,l}*x A 2 + S { 1 ,3, 1 } *y A 2 + S{l,l,3}*z A 2 + S{2,2,l}*x*y + S{2,l,2}*x*z +... 

S{ 1,2,2} *y*z + S{4,l,l}*x A 3 + S {3,2,1 }*x A 2*y + S{3,l,2}*x A 2*z + ... 

S{2,3,l}*x*y A 2 + S{l,3,2}*y A 2*z + S {2,1,3 }*x*z A 2 + S{l,2,3}*y*z A 2 +... 

S {2,2,2} *x*y*z + S{3,3,l}*x A 2*y A 2 + S{3,l,3}*x A 2*z A 2 + ... 

S{l,3,3}*y A 2*z A 2 + S {3,2,2} *x A 2*y*z + S {2,3,2} *x*y A 2*z + ... 

S {2,2,3 }*x*y*z A 2 + S{4,2,l}*x A 3*y + S{4,l,2}*x A 3*z + S{3,3,2}*x A 2*y A 2*z+... 

S {3,2,3 }*x A 2*y*z A 2 + S {2,3,3 }*x*y A 2*z A 2 + S{4,2,2}*x A 3*y*z + ... 

S {4,3,1 }*x A 3*y A 2 + S {4,1,3 }*x A 3*z A 2 + S{3,3,3}*x A 2*z A 2*y A 2 + - 

S {4,3,2} *x A 3*y A 2*z + S{4,2,3}*x A 3*y*z A 2 + S{4,3,3}*x A 3*y A 2*z A 2 ; 

LPV Representation 3 : Without using a For-loop, the equations were written without any particular order. 

LPV3 = S{3,l,l}*x A 2 + S{l,3,l}*y A 2 + S{l,l,3}*z A 2 + S{2,2,l}*x*y + S{2,l,2}*x*z +... 

S{ 1,2,2} *y*z + S{4,l,l}*x A 3 + S{3,2,l}*y*x A 2 + S{3,l,2}*z*x A 2 + ... 

S{2,3,l}*x*y A 2 + S{l,3,2}*z*y A 2 + S{2,l,3}*x*z A 2 +S{l,2,3}*y*z A 2+... 

S {2,2,2} *x*y*z + S{3,3,l}*(x A 2)*y A 2 + S{3,l,3}*(x A 2)*z A 2 + ... 

S{l,3,3}*(y A 2)*z A 2 + S {3,2,2} *(x A 2)*y*z + S {2,3,2} *x*y A 2*z + ... 

S {2,2,3 }*x*y*z A 2 + S {4,2,1 }*y*x A 3 + S{4,l,2}*z*x A 3 + S {3,3,2} *z*(x A 2)*y A 2 +•.. 

S {3,2,3 } *y*(x A 2)*z A 2 + S{2,3,3}*x*(y A 2)*z A 2 + S{4,2,2}*y*z*x A 3 + ... 

S {4,3,1 }*(x A 3)*y A 2 + S{4,l,3}*(x A 3)*z A 2 + S{3,3,3}*(x A 2)*z A 2*y A 2 + ... 

S {4,3,2} *z*(x A 3)*y A 2 + S{4,2,3}*y*(x A 3)*z A 2 + ... 

S {4,3,3 } *(y A 2)*(x A 3)*z A 2 ; 

Although the LPV1, LPV2, and LPV3 equations are exactly the same in terms of x, y, and z, the resulting LFT 
models (after reduction) were very different, as shown in Table 2. Based on the results for this example, LPV 
Representation 1 produced the LFT models of lowest dimension (and were therefore included in Table 1). Note that 
the ONERA tree-decomposition produced the same LFT model for each representation. The NASA method is also 
impervious to the equation format, because the LPV model does not need to be written in symbolic form - although 
re-defining the parameter order in the S matrix can have a minor impact on LFT model dimension for some 
problems. 
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Table 2. LFT Modeling Results Using ONERA and MuSyn Tools for Three Equivalent 
LPV Representations of the Generic Example 


\LPV Rep. 

LPV1 

LPV2 

LPV3 

Tool 

n A 

n x 

n y 

n z 

n A 

n x 

n y 

n z 

n A 

n x 

n y 

n z 

ONERA 

(“symtreed”) 













No Reduction 

110 

9 

24 

77 

110 

9 

24 

77 

110 

9 

24 

77 

1-D Reduction 

109 

9 

24 

76 

109 

9 

24 

76 

109 

9 

24 

76 

n-D Reduction 

46 

9 

24 

13 

46 

9 

24 

13 

46 

9 

24 

13 

ONERA 

(“lfrs”) 













No Reduction 

738 

318 

210 

210 

738 

318 

210 

210 

738 

318 

210 

210 

1-D Reduction 

445 

238 

148 

14 

448 

285 

135 

28 

471 

290 

135 

46 

n-D Reduction 

45 

9 

24 

12 

57 

9 

30 

18 

96 

45 

27 

24 

MuSyn 

(“ureal”) 













No Reduction 

738 

318 

210 

210 

738 

318 

210 

210 

738 

318 

210 

210 

Basic 

579 

268 

173 

138 

632 

268 

182 

182 

608 

270 

169 

169 

Full 

45 

9 

24 

12 

57 

9 

30 

18 

96 

45 

27 

24 


5. Conclusion 

Methods and software tools for developing linear fractional transformation (LFT) models for uncertain systems 
were considered in this paper, as a precursor to applying formal robustness analysis methods to control upset 
prevention and recovery systems as part of a validation process being developed for potential use in their ultimate 
certification. Such systems (developed for failure detection, identification, and reconfiguration, as well as upset 
recovery) need to be evaluated over broad regions of the flight envelope and under extreme flight conditions, and 
should include various sources of uncertainty. Flowever, formulation of LFT models for representing system 
uncertainty can be very difficult for complex parameter-dependent systems. A numerical matrix-based LFT 
modeling method and preliminary software tool were presented and evaluated in this paper for several example 
problems and in comparison to other available methods and software tools. The examples included an F-16 aircraft 
at an extreme flight condition (i.e., near a stall bifurcation), a missile model problem, and a generic example 
designed to be a difficult modeling problem (particularly relative to cross-product terms). The numerical modeling 
method and preliminary software tool presented in this paper compared favorably for each of the example problems 
relative to the other methods considered. The matrix-based modeling approach therefore appears to be promising. 
Several areas for further refinement of the preliminary software tool were also discussed. Further research will 
focus on these refinements, as well as applying the tool to robustness analysis studies for control upset prevention 
and recovery technologies. These studies will provide risk mitigation for high-risk flight testing under extreme 
and/or loss-of-control flight regimes, aircraft failure and damage, and other adverse or upset conditions. These kinds 
of high-risk tests will be performed at NASA Langley using a dynamically scaled transport aircraft model, as part of 
the Airborne Subscale Transport Aircraft Research (AirSTAR) Testbed. Open-loop tests will be performed to 
validate and further investigate vehicle dynamics under extreme/upset conditions, and closed-loop tests will be 
performed for failure accommodation, upset recovery, and damage mitigation. 

A possible advantage of the numerical LFT modeling method presented in this paper is its potential future use as 
part of an online robustness analysis tool for risk mitigation during high-risk flight tests or for onboard aircraft 
applications. An on-line robustness analysis tool is currently being developed to provide risk mitigation during 
flight tests involving the AirSTAR Testbed, and a future extension could possibly utilize an online uncertainty 
modeling capability to update the system model being used for analysis. Onboard modeling and robustness analysis 
methods for future transport aircraft applications may also be feasible. 
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Appendix 


The LPV model for the F-16 example presented in Section 3 is given as follows. 

A(s[,.? 2 ) = AAOO + AAlOs t + AAOls-, + AAl Is,,?, + AA20s 2 
+AAQ2s 2 + AA2\sls 2 + AA\2s l s 2 + AA30s^ + AA03s 2 

B (.S', ,.S' 2 ) = BBOO + BB\0s l +BB01S-, + .5.51 + BB20s* 

+BB02s; + BB2ls 2 s 2 + BB12 s/ 2 + BB30S* + BB03s 2 2 

C(s t ,s 2 ) = CCOO + CCIO^ + CC0b 2 + CC1 Is,.?, + CC20s 2 
+CC02s; + CC2ls 2 s 2 + CCl2s t s 2 2 + CC30^ 3 + CC03s 2 

D(s t ,s 2 ) = DD00 + DDlOs l +DD0ls 2 + .DZ)1 +DD20s 2 

+DD02s 2 + DD2ls 2 s 1 +Z)Z)12s , 1 s' 3 + DD30s 2 +DD03sl 


The nonzero coefficient matrices are given below. 


0 0 0 

0 0 0 

0 0 0 

0 0 0 

AAOO =o 0 0 

0 0 0 

0 -4.1662 0 

0 0.0003 0 

0.0206 0 0 


1.0000 0 1.1702 

0 1.0000 0 
0 0 1.5393 

-0.0042 0 -0.0849 

0 -0.0152 0 

-0.0001 0 -0.0131 

0 -10.2531 0 

0 0.9241 0 

0.7607 0 -0.6589 


0 0 0 

0 0 0 

0 0 0 

0 0 0.0157 

0 0.0127 0 

0 0 -0.0668 

0.2358 -9.4733 0 

0 0.1824 0 

0 0 -0.0822 


0 

0 

0 

0 

0 

-0.0251 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-0.0191 

0 

0 

0 

0 

0 

0 

-0.0002 

0 

0.0088 

0 

0 

0.0016 

0 

0 

0 

0 

-0.0023 

0 

0 

-0.0006 

0 

0 

0 

0 

0.0003 

0 

0.0003 

0 

0 

0.0016 

0 

0 

0 

0 

-0.1232 

0 

-0.0003 

-2.3492 

0 

0 

0 

0 

0 

-0.0013 

0 

0 

-0.0151 

0 

0.0003 

0 

0 

-0.0056 

0 

-0.0063 

0 

0 

0.0002 
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. 4.402 = 10 4 


. 4.403 = 10 : 


. 4.410 


AA11 = 



0 

0 

0 

0 

0 

- 0.0779 

0 

0 

0 



0 

0 

0 

0 

0 

0 

0 

0 

0 



0 

0 

0 

0 

0 

- 0.0506 

0 

0 

0 



0 

0 

0 

- 0.0018 

0 

0.0159 

0 

0 

- 0.0151 


0 

0 

0 

0 

0.0009 

0 

0 

0.0019 

0 



0 

0 

0 

0.0002 

0 

0.0006 

0 

0 

0.0034 


0 

0.0003 

0 

0 

0.6992 

0 

0.0096 

- 3.9688 

0 



0 

0.0001 

0 

0 

- 0.0019 

0 

- 0.0002 

- 0.0543 

0 


L 0.0004 

0 

0 

- 0.0300 

0 

- 0.0272 

0 

0 

-O.OOigJ 



0 

0 


0 0 

0 

- 0.3310 

0 

0 

0 

“ 



0 

0 


0 0 

0 

0 

0 

0 

0 




0 

0 


0 0 

0 

- 0.1872 

0 

0 

0 




0 

0 


0 -0.0110 0 

0.0055 

0 

0 

- 0.1278 



0 

0 


0 0 

0.0256 

0 

0 

- 0.0059 

0 




0 

0 


0 - 0.0016 0 

0.0003 

0 

0 

0.0083 



0 

0.0025 

0 0 

5.8716 

0 

0.0869 

- 4.0137 

0 




0 

0.0011 

0 0 

0.0115 

0 

- 0.0006 

- 0.1937 

0 




- 0.0008 0 


0 - 0.1537 0 

- 0.1127 

0 

0 

- 0 . 0121 . 


0 

0 

0 

0 

0 

0.2958 

0 

0 

0 



0 

0 

0 

0 

0 

0 

0 

0 

0 



0 

0 

0 

0 

0 

0.2249 

0 

0 

0 



0 

0 

0 

0.0030 

0 

- 0.1147 

0 

0 - 0.0212 


= 


0 

0 

0 

0 

0.0304 

0 

0 0.0082 

0 



0 

0 

0 

- 0.0035 

0 

- 0.0039 

0 

0 - 0.0218 



0 

0 

0 

0 

1.5240 

0 0.0016 30.8244 

0 



0 

0 

0 

0 

0.0173 

0 0.0005 0.2009 

0 



- 0.0032 0 

0 

0.0722 

0 

0.0815 

0 

0 - 0.0007 


0 


0 

0 

0 

0 

18.3487 

0 

0 

0 

-1 

0 


0 

0 

0 

0 

0 

0 

0 

0 


0 


0 

0 

0 

0 

11.9117 

0 

0 

0 


0 


0 

0 

0.4309 

0 

- 3.1205 

0 

0 

4.1293 

0 


0 

0 

0 

- 0.4923 

0 

0 

- 0.5664 

0 


0 


0 

0 

- 0.0330 

0 

- 0.1094 

0 

0 

- 0.6622 

0 

- 

- 0.0196 

0 

0 

- 194.1753 

0 

- 2.4965 

740.3101 

0 


0 

- 

- 0.0129 

0 

0 

0.3508 

0 

0.0413 

12.3560 

0 


0.0864 

0 

0 

7.2006 

0 

6.4003 

0 

0 

0.5259 I 
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AA12 


. 4.420 


AA21 = 


. 4.430 


0 

0 

0 

0 

0 

0.1169 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0.0661 

0 

0 

0 

0.0040 

0 

- 0.0008 

0 

0 

0 

0 

- 0.0092 

0 

0 

0 

0 

0.0006 

0 

- 0.0001 

0 

- 0.0008 

0 

0 

- 2.1076 

0 

0 

- 0.0004 

0 

0 

- 0.0046 

0 

0.0004 

0 

0 

0.0548 

0 

0.0396 


0 0 0 

0 0 0 

0 0 0 

0 0 0.0460 

0 0.0029 0 

0 0 - 0.0029 

- 0.0317 1.2821 0 

0.0002 0.0675 0 

0 0 0.0044 


0 

0 

0 

0 

0 

- 0.1080 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

- 0.0701 

0 

0 

0 

0 

0 

0 

- 0.0026 

0 

0.0137 

0 

0 

- 0.0281 

0 

0 

0 

0 

0.0049 

0 

0 

0.0042 

0 

0 

0 

0 

0 

0 

0.0005 

0 

0 

0.0028 

0 

0.0001 

0 

0 

1.3392 

0 

0.0161 

- 2.9008 

0 

0 

0.0001 

0 

0 

- 0.0014 

0 

- 0.0002 

- 0.0687 

0 

0.0004 

0 

0 

- 0.0430 

0 

- 0.0373 

0 

0 

- 0.0035 


0 


0 0 

0 

0 

- 13.7634 

0 

0 

0 

0 


0 0 

0 

0 

0 

0 

0 

0 

0 


0 0 

0 

0 

- 7.7846 

0 

0 

0 

0 


0 0 

- 0.4737 

0 

0.0278 

0 

0 

- 5.4545 

0 


0 0 

0 

1.0739 

0 

0 

- 0.4604 

0 

0 


0 0 

- 0.0765 

0 

0.0050 

0 

0 

0.3573 

0 


0.0782 0 

0 249.1658 

0 

3.8505 

- 159.6116 0 

0 


0.0452 0 

0 

0.5933 

0 

- 0.0243 

- 7.8976 

0 

- 0.0481 

0 0 

- 6.5172 

0 

- 4.6449 

0 

0 

- 0.5303 


0 

0 

0 0 

0 

0.0540 

0 

0 

0 



0 

0 

0 0 

0 

0 

0 

0 

0 



0 

0 

0 0 

0 

0.0306 

0 

0 

0 



0 

0 

0 0.0019 

0 

- 0.0002 

0 

0 

0.0213 


= 

0 

0 

0 0 

- 0.0040 

0 

0 

0.0024 

0 



0 

0 

0 0.0003 

0 

0 

0 

0 

- 0.0016 



0 

- 0.0003 

0 0 

- 0.9675 

0 

- 0.0155 

0.7910 

0 



0 

- 0.0002 

0 0 

- 0.0025 

0 

0.0001 

0.0312 

0 



0.0002 0 

0 0.0258 

0 

0.0182 

0 

0 

0.0021 
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BB 00 = 


BB 01 


BB 02 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-0.0167 

-0.3934 

0.3934 

0.0045 

0 

0 

-0.0111 

-0.0111 

0 

0 

0.0009 

-0.0061 

0.0061 

-0.0014 

0 

0 

-2.3330 

-2.3330 

0 

0.0010 

0 

-0.0038 

-0.0038 

0 

0 

0.0044 

0 

0 

0.0121 

0 


0 0 

0 0 

0 0 

-0.4317 0 

0 -0.8317 

0.1155 0 

0 20.8448 

0 -0.3167 

0 0 


0 0 0 

0 0 0 

0 0 0 

0 0.1210 0 

-0.8317 0 0 

0 -0.1519 0 

20.8448 0 0.0116 

-0.3167 0 0.0001 

0 0 0 


0 0 0 

0 0 0 

0 0 0 

-50.6223 -214.8311 214.8311 

0 -12.7580 -12.7580 

1.5743 -3.3433 3.3433 

0 -645.2558 -645.2558 

0 -8.5432 -8.5432 

1.2085 0 0 


0 0 

0 0 

0 0 

8.0884 0 

0 0 

-6.7643 0 

0 0.4809 

0 0.0063 

3.2993 0 
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BB 03 = 


BB 10 


BBU 



0 

0 

0 

0 

0 


0 

0 

0 

0 

0 


0 

0 

0 

0 

0 


-0.3677 

-1.7900 

1.7900 

0.0527 

0 

4 







0 

-0.2239 

-0.2239 

0 

0 


-0.0019 

-0.0279 

0.0279 

-0.0371 

0 


0 

-3.9126 

-3.9126 

0 

0.0019 


0 

-0.0678 

-0.0678 

0 

0 


. 0.0101 

0 

0 

0.0275 

0 

r 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0.0054 

-0.0060 

0.0060 

-0.0015 

0 


0 

0.0106 

0.0106 

0 

0 

- 

0.0015 

-0.0001 

0.0001 

0.0020 

0 


0 

-0.3065 

-0.3065 

0 

-0.0002 


0 

0.0041 

0.0041 

0 

-0 

- 

0 

0 

0 

0.0001 

0 

r 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

1.2691 

5.7140 

-5.7140 

-0.1971 

0 


0 

0.2341 

0.2341 

0 

0 

- 

0.0272 

0.0889 

-0.0889 

0.1578 

0 


0 

20.0750 

20.0750 

0 

-0.0113 


0 

0.1866 

0.1866 

0 

-0.0002 


0.0321 

0 

0 

-0.0878 

0 
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0 0 0 0 0 
0 0 0 0 0 
0 0 0 0 0 


0.1329 0.6475 

BB12 = 10 3 o 0.0865 
0.0008 0.0101 
0 1.3147 

0 0.0253 

.-0.0037 0 

0 0 

0 0 

0 0 

-0.0080 -0.0383 
BB20 = o -0.0009 
0.0001 -0.0006 

0 -0.1536 

0 -0.0010 

. 0.0002 0 

0 0 

0 0 

0 0 

-1.5956 - 7.7385 

BB21 = o -1.1246 
-0.0097 -0.1204 

0 -13.9688 

0 -0.3181 

0.0438 0 


-0.6475 -0.0190 0 

0.0865 0 0 

-0.0101 0.0133 0 

1.3147 0 -0.0007 

0.0253 0 0 

0 - 0.0100 0 

0 0 0 

0 0 0 

0 0 0 

0.0383 0.0012 0 

-0.0009 0 0 

0.0006 -0.0009 0 

-0.1536 0 0.0001 


- 0.0010 0 0 


0 

0.0006 

o J 

0 

0 

0 

0 

0 

0 

0 

0 

0 

7.7385 

0.2281 

0 

-1.1246 

0 

0 

0.1204 

-0.1599 

0 

-13.9688 

0 

0.0078 

-0.3181 

0 

0.0002 

0 

0.1195 

0 
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BB20 = 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0.0064 

0.0305 

- 0.0305 

- 0.0009 

0 

0 

0.0049 

0.0049 

0 

0 

0 

0.0005 

- 0.0005 

0.0006 

0 

0 

0.0453 

0.0453 

0 

0 

0 

0.0013 

0.0013 

0 

0 

- 0.0002 

0 

0 

- 0.0005 

0 


0 0 0 0 0 0 0 
0 0 0 0 0 0 0 
1 0 0 0 0 0 0 
0 1 0 0 0 0 0 
0 0 1 0 0 0 0 
0 0 0 1 0 0 0 
0 0 0 0 1 0 0 
0 0 0 0 0 1 0 
0 0 0 0 0 0 1 
0 0 0 0 0 1 0 


ccoo = 


0 0 
0 0 
0 0 
0 0 


0 0 
0 -1 
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